**** Synthetic control placebo testing 

matrix deterrence = (0,0)

*** Load in data, loop over different controls
cd /disk/agedisk3/medicare.work/poterba-DUA52260/jetson-dua52260/kyphon/synthcontrolsshort/100pct
log using synth_placebo.log, replace
use growthgroupspanel.dta, clear 
drop if group == 999
levelsof group, local(placebos)

foreach p of local placebos{


	cd /disk/agedisk3/medicare.work/poterba-DUA52260/jetson-dua52260/kyphon/synthcontrolsshort/100pct/placebos

	*** Grab full data and remove treated unit
	use /disk/agedisk3/medicare.work/poterba-DUA52260/jetson-dua52260/kyphon/synthcontrolsshort/100pct/growthgroupspanel.dta, clear 
	drop if group == 999 

	***********************************
	* Swap in placebos
	disp "Placebo group `p'"
	replace group = 999 if group == `p'
	save growthgroupspanel_temp.dta, replace 


	*********************************************************************
	*** Regress Placebo on lags to find best fit  

	levelsof group, local(levels)
	keep group dt pmt_amt 
	reshape wide pmt_amt, i(dt) j(group)
	tsset dt
	* 3 columns: group, lag period, distance
	mata: distanceLead= J(1,3,0)
	mata: distanceLag = J(1,3,0)


	set matsize 11000

	** Only want to keep kyphoplasty pre-period 
	local treatment_period =  mofd(mdy(12, 12, 2005))
	replace pmt_amt999 = . if dt >= `treatment_period'

	** Find max lead so that there are 5 years of post-period data 
	summarize dt 
	local maxlead = r(max) - `treatment_period' + 1 - 60  
		** Max - treatment  + 1 = length of end window 
		** Take 60 off and that's the max you can shift before running out of data

	** Find max lag so that there are 3 years of pre-period data 
	summarize dt 
	local maxlag = `treatment_period' - r(min) - 36
		** Treatment period - min  = length of end window 

	gen diff = 0 

	*** Loop, storing leads lags and fit
	foreach g of local levels {
		disp `g'
		* Skip Kypho group
		if `g' == 999{
				continue
			}

	foreach t of numlist 0/`maxlag'{
		*disp `t'
		qui replace diff =  (pmt_amt999 - L`t'.pmt_amt`g')^2
		qui summarize diff 
		local dist = r(sum)/r(N)
		mata: distanceLag=  distanceLag \ `g', `t', `dist'
		}

	foreach t of numlist 0/`maxlead'{
		* disp `t'
		qui replace diff = (pmt_amt999 - F`t'.pmt_amt`g')^2
		qui summarize diff 
		local dist = r(sum)/r(N)
		mata: distanceLead=  distanceLead \ `g', `t', `dist'
		}
	}

	mata: st_matrix("distanceLead", distanceLead)
	mata: st_matrix("distanceLag", distanceLag)

	svmat distanceLag
	keep distanceLag*
	svmat distanceLead

	rename distanceLag1 laggroup
	rename distanceLag2 lagperiod
	rename distanceLag3 lagdistance
	rename distanceLead1 leadgroup
	rename distanceLead2 leadperiod
	rename distanceLead3 leaddistance

	drop if laggroup == 0
	* Don't save, immediately use 
	********************************************************************************
	*** Find best fit per group

	preserve
	keep laggroup lagperiod lagdistance
	drop if mi(laggroup)
	bysort laggroup: egen minlagdistance = min(lagdistance)
	keep if lagdistance == minlagdistance
	drop minlagdistance
	rename laggroup group
	sort group lagperiod 
	collapse (first) lagperiod (first) lagdistance, by(group)
	save group_lag_mindistance_temp.dta, replace

	restore
	keep leadgroup leadperiod leaddistance
	drop if mi(leadgroup)
	bysort leadgroup: egen minleaddistance = min(leaddistance)

	keep if leaddistance == minleaddistance
	drop minleaddistance
	rename leadgroup group
	sort group leadperiod 
	collapse (first) leadperiod (first) leaddistance, by(group)
	save group_lead_mindistance_temp.dta, replace

	clear
	use group_lag_mindistance_temp.dta
	merge 1:1 group using group_lead_mindistance_temp.dta

	drop _merge
	bysort group: gen mindistance = min(leaddistance, lagdistance)

	replace lagdistance = . if lagdistance!= mindistance
	replace leaddistance = . if leaddistance!= mindistance

	replace lagperiod = . if lagdistance == .
	replace leadperiod = . if leaddistance == .

	replace lagperiod = -1*lagperiod
	gen period = min(lagperiod, leadperiod)
		* Min because missing is infty 

	keep group period mindistance
	drop if mi(group)
	save groupleadlags_temp.dta, replace

	************************************
	*** Run Synthetic controls on fit lead/lags

	cd /disk/agedisk3/medicare.work/poterba-DUA52260/jetson-dua52260/kyphon/synthcontrolsshort/100pct/placebos
	use growthgroupspanel_temp.dta, clear

	*** Merge leads and lags
	merge m:1 group using groupleadlags_temp.dta
	drop _merge mindistance


	** Store group levels in macro
	levelsof group, local(levels)


	*** Replace generate with fitted lead/lags
	sort group dt
	gen out = .


	foreach g of local levels{
		disp `g'
		* Treated group gets no time shift
		if `g' == 999{
			replace out = pmt_amt if group == `g'
			continue
			}

		** Grab the stored shifter 
		qui summarize period if group == `g'
		local shift = r(min)

		** If positive, create lead 
		if `shift' > 0{
			replace out = F`shift'.pmt_amt if group == `g'
		}
		** If negative, create lag 
		if `shift' < 0{
			local shiftvalue = abs(`shift')
			replace out = L`shiftvalue'.pmt_amt if group == `g'
		}

		if `shift' == 0{
			replace out = pmt_amt if group == `g'
		}
	}


	** Cut missing data from front of time series
	** Which is limited by the highest lag 
	summarize period
	drop if dt < mofd(mdy(1,1,1999)) + abs(min(r(min), 0))

	*** Cut missing data from end of time series
	*** Which is limited by furthest lead
	*** Data end 9/1/2015 because of second DRG change
	summarize period
	drop if dt > mofd(mdy(9, 1, 2015)) - abs(max(0,r(max)))


	*** Filing date
	disp mofd(mdy(12, 12, 2005))
		* 551 



	*** We're using filing date 
	format dt %tm_Mon_YY

	**** Run Synthetic Controls
	synth out out, trunit(999) trperiod(551) figure

	*** Store fitted values
	matrix fit = e(Y_synthetic)
	keep if group == 999
	svmat fit 
	save synth_fit_temp.dta, replace

	********************************************************************************
	*** Compute deterrence
	cd /disk/agedisk3/medicare.work/poterba-DUA52260/jetson-dua52260/kyphon/synthcontrolsshort/100pct/placebos
	use growthgroupspanel_temp.dta, clear
	keep if group == 999 

	merge 1:1 dt using synth_fit_temp.dta


	local treatment_period =  mofd(mdy(12, 12, 2005))
	disp `treatment_period'

	** Produce difference in series
	gen deterrence = fit1 - out if dt >= `treatment_period'


	** Time discount the difference at 10% 
	local r = 1.1^(1/12)

	gen deterrence_disc = deterrence*1/(`r'^(dt-`treatment_period'))
	replace deterrence_disc = . if dt-`treatment_period' >= 60
	count if !mi(deterrence_disc)
	assert(r(N)==60)

	*** Compute total
	collapse (sum) deterrence_disc
	summarize deterrence_disc
	local deterrence = r(min)
	
	matrix deterrence = (deterrence\ `p', `deterrence')

	rm synth_fit_temp.dta 
	rm groupleadlags_temp.dta 
	rm group_lead_mindistance_temp.dta 
	rm group_lag_mindistance_temp.dta 
	rm growthgroupspanel_temp.dta 

}

clear 
svmat deterrence
rename deterrence1 group 
rename deterrence2 deterrence
sort deterrence
save deterrence_dist.dta, replace

*** Deterrence computation
cd /disk/agedisk3/medicare.work/poterba-DUA52260/jetson-dua52260/kyphon/synthcontrolsshort/100pct/placebos
use deterrence_dist.dta, clear
drop if group == 0
count 
	* 30
count if deterrence > 5.389*10^8
	* 4

log close 

